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Abstract. This paper concerns the reconstruction of a scalar diffusion coef- 
ficient cr(x) from redundant functionals of the form Hi(x) = <t 2q (x)\ Vu; \ 2 (%) 
where a S K and iz; is a solution of the elliptic problem V • crVu; = for 
1 < i < J. The case a = | is used to model measurements obtained from 
modulating a domain of interest by ultrasound and finds applications in ul- 
trasound modulated electrical impedance tomography (UMEIT), ultrasound 
modulated optical tomography (UMOT) as well as impedance acoustic com- 
puterized tomography (ImpACT). The case a = 1 finds applications in Mag- 
netic Resonance Electrical Impedance Tomography (MREIT). 

We present two explicit reconstruction procedures of a for appropriate 
choices of / and of traces of Ui at the boundary of a domain of interest. The 
first procedure involves the solution of an over-determined system of ordinary 
differential equations and generalizes to the multi-dimensional case and to (al- 
most) arbitrary values of a the results obtained in two and three dimensions 
in [10] and [5], respectively, in the case a = ~. The second procedure consists 
of solving a system of linear elliptic equations, which we can prove admits a 
unique solution in specific situations. 

1. Introduction. Medical imaging modalities aim to combine high resolution with 
high contrast between healthy and unhealthy tissues. Optical Tomography and 
Electrical Impedance Tomography display such high contrasts but often suffer from 
poor resolution. Ultrasound Tomography and Magnetic Resonance Imaging are 
high resolution modalities that sometimes suffer from low contrast. The ultrasound 
modulation of electrical or optical properties of tissues and the combination of simul- 
taneous electrical and magnetic resonance measurements both offer the possibility 
to combine high resolution with high contrast. For the acquisition of ultrasound- 
modulated measurements, we refer the reader to, e.g., [1, 5, 10, 12, 14] for works 
in the mathematical literature. For the acquisition of internal information on elec- 
trical conductivities by magnetic resonance imaging, we refer the reader to, e.g., 
[13, 17, 18, 19]. 
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Mathematically, we aim to reconstruct a scalar diffusion coefficient a in an 
elliptic equation from knowledge of internal information of the form Hij(x) = 
a 2a (x)\7ui(x) ■ S7uj(x) for a G R and 1 < i,j < m, where and Uj are solu- 
tions of the elliptic problem with different boundary conditions; see (1) and (2) 
below. Coupling impedance (or diffusion) with acoustic waves or magnetic reso- 
nance correspond to the cases a = \ and a = 1, respectively. Such information can 
be obtained from functionals of the form a 2a (x)\ s \7ui\ 2 (x) by standard polarization 
(expressions of the form Aab = (a + 6) 2 — (a — b) 2 ). 

This problem was first solved in the two dimensional setting in [10] in the case 
m = 2 and a = i. The three dimensional setting was addressed in [5] with m = 4 
and a = i. In these papers, the elliptic equation is recast as a system of equations 
for quantities of the form Si = a a X7ui using the elliptic equation and the fact that 
Vui is curl free. This strategy allows one to eliminate a from the system of equations 
and solve for the vectors Si . The stable reconstruction of a is then straightforward. 
The case a = | in the setting of non- redundant measurements, i.e., with m = 1 and 
measurements of the form H = ct|Vu| 2 is considered in [4]. It is shown in that paper 
that the stable reconstruction of a may not be possible from such non-redundant 
measurements. This justifies the analysis of redundant measurements. 

The objectives of this paper are twofold. We first generalize the reconstruction of 
a to the case of arbitrary space dimension n and almost arbitrary a € R. Assuming 
that the vectors Si form a frame, we obtain a system of equations for the vectors 
Si that involves Hij = Si ■ Sj but no longer a. The resulting system of equations 
may be seen as an overdetermined nonlinear system of equations. By appropriately 
choosing the boundary conditions used to construct the internal functionals Hij{x), 
we obtain a global uniqueness and stability result for the reconstruction of the 
scalar quantity a(x). Although several portions of the algorithm generalize to the 
reconstruction of anisotropic diffusion tensors, we restrict ourselves to the scalar 
case in this paper. We also describe and investigate the compatibility conditions 
associated with such a redundant system. 

The second objective of the paper is to present a system of elliptic equations 
for the solutions m with constitutive parameters that depend on the measurements 
Hij but not on the unknown diffusion coefficient a. We show that the system is 
uniquely solvable when a Fredholm alternative holds. We obtain existence and 
uniqueness results for the proposed system for all but a discrete number of values 
of the dimension n and the coefficient a € R. 

Both algorithms require boundary conditions for the elliptic solutions that ensure 
that n of the vectors Si form a frame in R" at each point of the domain of interest. 
Whereas such a condition is easy to meet in two space dimensions, in dimensions 
three and higher, the only available technique that guarantees such an independence 
is based on using complex geometrical optics (CGO) solutions. We generalize here 
the CGO construction of [5] to the multi-dimensional setting and for almost all 
values of a. 

The inverse diffusion problems with internal functionals considered here are ex- 
amples of hybrid inverse problem where two imaging modalities are combined to 
provide both high resolution and high contrast. For recent works on the mathemat- 
ics of hybrid inverse problems and their many applications in medical imaging, we 
refer the reader to the articles in the book [20] and to the recent review paper [3]. 

The rest of the paper is structured as follows. Section 2 presents the main results 
of the paper on the stable reconstruction of a from available internal functionals. 
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The elimination of a from the system of equations for the vectors Si and the cor- 
responding differential calculus is explained in section 3. The redundant system 
of equations for the vectors Si is addressed in section 4 while the system of linear 
equations for the solutions Ui is given in section 5. Finally section 6 presents further 
reconstruction algorithms in the two dimensional case and analyzes the compatibil- 
ity conditions satisfied by the redundant data and their potential use. 

2. Statement of the main results. Let X be an open convex bounded domain 
of M. n with n > 2. In the following, we address the reconstruction of the scalar 
conductivity (or diffusion) coefficient a in the equation 

V-{aVui)=0 (X), Ui\ ax =9i, l<i<m, (1) 

where m > n, from knowledge of the interior functionals 

H i:j (x) = a{x) 2a Vui(x) ■ Vu.j(x), l<j<i<m, (2) 

where aglis fixed and such that (n — 2)a + 1^0. The derivation of the internal 
functionals (2) in the case a = | is detailed in [5, 14] as examples of synthesized 
focusing, in [1] in a setup of temporal, physical, focusing, and in [12] by considering 
thermoelastic effects. The case a = 1 with m = 1 related to MREIT and CDII 
(Current Density Impedance Imaging) is addressed in [13, 17, 18, 19]. 

Following a similar approach to [5, 10], we first perform the change of unknown 
functions Si = <r a Vui for every i and define 

F(x) := V\oga(x). (3) 

We also equip X C R™ with its Euclidean metric = in the canonical basis 
(ei, . . . , e„). For a given vector field V = V l ei defined on X, we define the cor- 
responding one-form V := V l dx l (i.e., by means of the flat operator). With this 
notation, we obtain that the vector fields Sj satisfy the system of equations 

V ■ Sj = -(l-a)F ■ Sj, (4) 

dS) = aF b A Sj, 1 < j < m, (5) 

where A and d denote the usual exterior product and exterior derivative, respec- 
tively. The first equation stems directly from (1) whereas the second one states 
that the one-form o-~ a S b j = duj is exact, therefore closed, and hence d(a~ a S b j) = 0. 
When n = 2, 3, equation (5) is recast as: 

n = 2 : V 1 - ■ Sj - aJF -Sj = 0, n = 3 : curl Sj - aF x Sj = 0, 

where in dimension n = 2, we define J := [° ~q] and V- 1 := JV, and in dimension 
n = 3, x denotes the standard cross-product. The available information becomes 
Hij(x) = Si(x) ■ Sj(x). 

A crucial hypothesis for our reconstruction procedure is that the m gradients 
have maximal rank in M™ at every point x € X. This hypothesis can be for- 
malized by the somewhat stronger statement: there exists a finite open covering 
O = {Qk}i<k<N of X (i.e. X C U^ =1 flk), an indexing function r : [1,N] 3 i i-> 
t(i) — (r(i)i, . . . ,t(z)„) £ [l,m] n and a positive constant Cq such that 

^min^^mf det(S , T ( i)l (x), . . . , S T ^) n (x)) > c > 0. (6) 

This assumption is equivalent to imposing the following condition on the data 

man inf det H r{l) (x) > cl > 0, (7) 

l<i<Nx£Qi 
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where H T ^ stands for the nxn matrix of elements = 5V(i) fc • Sr(i)i ■ While one 

can always find illuminations such that (6) holds in two dimensions with m = n = 2 
and O = {X} (the most prcferrablc case) by virtue of [2, Theorem 4], higher 
dimensions can be dealt with using complex geometrical optics solutions provided 
that a has enough regularity, as the following lemma shows. 

Lemma 2.1. Let n > 3 and a G H% +3+£ (X) for some e > be bounded from below 
by a positive constant. Then 

(i) for n even, there exists a non-empty open set G of illuminations {g\, . . . , g n } 
such that for any g € G, the condition (7) holds with O — {X} for some 
constant cq > 0. 

(ii) For n odd, there exists a non-empty set G of illuminations {g\, .., g n +i} such 
that for any g G G there exists an open cover of X of the form {Q2i-1) ^2i}i<i<./v 
and a constant Co > such that 

inf det(Si, . . . , S n _x,£iS n ) > c and inf det(5i, . . . , S n _i, e^n+i) > 0) = 

(8) 

for 1 < i < N and with £i,£i — ±1. 

The first step toward an inversion is to express the source term F in terms of a 
local frame: 



Lemma 2.2. Let Q C X open where, up to renumbering solutions, we have 
inf dct(S'i(a;), . . . , S n (x)) > cq > 0. 

Then at every point x € and denoting H(x) := {Si(x) ■ Sj(x)}i<ij< n , D{x) = 
y det H(x), the vector field F(x) = Vlogcr(x) is given by the following formulas 

(9) 



F = % E WDHV) ■ Si)Si = c F (V log D + J2 ( Vffy ' Si)Si) , 



c F := ((n-2)a + l)- 1 . 
where denotes the element (i,j) of the matrix H^ 1 . 

Formula (9) was first proved in [5] in the two- and three-dimensional cases with 
a = 2 an d is here proved for general n and a e E such that (n — 2)a + 1^0. 
This formula gives us a way to reconstruct F locally from n linearly independent 
solutions. Assuming condition (7), one can then reconstruct F globally over X . 

From lemma 2.2, one can follow two directions to reconstruct the conductivity, 
which we now describe in more detail in the next two paragraphs. 
The ODE-bascd reconstruction procedure. The first approach consists in plugging 
equation (9) back into the system (4)-(5) and obtain a closed system for the vectors 
Sj. We then show that the resulting system leads to a gradient system, which 
can then be solved for the vectors Sj by ODE integration. Once the vectors Sj are 
reconstructed, one recovers a from the knowledge of its value at a given point and the 
fact that Vlogcr is now known by equation (9). This approach is a generalization 
of the results of [5] to higher-dimensional settings and general a G K such that 
(n — 2)q + 1^0, and leads to well-posed reconstructions as stated in the following: 

Theorem 2.3 (Global uniqueness and stability, ODE-based reconstruction proce- 
dure). Let X C > 2 be an open convex bounded domain, and let two sets of 
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m>n solutions of (1) generate measurements (H,H') whose components belong to 
W ,00 (X), and who jointly satisfy condition (7) with the same triple (0,t,cq). Let 
also x e fij C X and o-(x ),a' (x ) and {<SV(i )i (x ), (xo)}i<i<n fee gwen. 

Le£ cr and cr' fee £/ie conductivities corresponding to the measurements H and H' , 
respectively. Then we have the stability estimate: 

|| logo- - logt/HwrLoopt) <C(e + \\H- H'\\ W i,oo {x) ) , (10) 

where Eq is the error committed at the point Xq : 

n 

£ := |logcr(x ) - logCT / (x )| +^2\\S T ( io ) i (x ) - S' T{lo)z (x )\\. 

i=l 

The solution for the vectors Sj and then for logcr requires the solution of full 
gradient equations of the form Vw = f{u) 1 where u stands for either unknown. 
These overdetermined PDEs require compatibility conditions on / if we wish to 
ensure that their solution does not depend on the path of integration. Theorem 2.3 
shows that the reconstruction is unique and stable with respect to the data once a 
fixed family of integration curves is chosen. The compatibility conditions addressed 
in section 6 are shown to depend quadratically on the unknown frame. It is therefore 
difficult to enforce them while solving for the frame S. Nonetheless, depending on 
the value of a, they may lead to algebraic (i.e. pointwise) reconstructions of all or 
part of the unknown frame, and may also provide further conditions on the data 
Hij. Such analyses are carried out in section 6. 

Remark 1. Solving a system of equations for the unknown Si may not be efficient 
numerically. Let S be the matrix whose columns are the n linearly independent 
vectors Sj at a given x. Then S T S = H is known. By the Gram-Schmidt (GS) 
orthonormalization procedure or by setting R — SH~^ , we can write an equation 
for an oriented orthonormal frame R; see section 4.3 below. This approach requires 
that we reconstruct n(n — l)/2 = dim SO n (W) scalar functions instead of the nx m 
components of the vector fields {Sj}. The only additional constraint is that the 
transition matrix from S to R satisfies a certain stability property with respect to 
the data H, see Section 4.3 for details. 

Remark 2. The case a = 0, corresponding to information of the form Hij(x) = 
\7ui(x) ■ Vitj(x) = Si(x) ■ Sj(x) simplifies in the sense that the elimination of F is 
not necessary. Indeed, we show in the next section that knowledge of Si ■ Sj and 
the constraints dS^ = for 1 < j < m uniquely determine the vectors Sj provided 
they are known at one point xq. Once V«i is known, the reconstruction of a may 
proceed from using (9). Note that, alternatively, the equation (1) may be seen as a 
transport equation for cr 1_Q when a ^ 1 once the vector field er a Vit is known. The 
stability properties of such a reconstruction are established in [6, 7]. 

The elliptic-based reconstruction procedure. The second approach is novel and con- 
sists in injecting equation (9) back into the initial conductivity equations and obtain 
a strongly coupled elliptic system of the form 

Au t + c F Wij ■ Vuj = 0, Ui\ dx — g%, l<i<m, (11) 

where the vector fields Wij are known from the data and where the illuminations gi 
were prescribed in the first place. Here and below, we use the Einstein convention of 
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summation over repeated indices. The vector fields Wij satisfy stability conditions 
of the form 

||W||oo := max ||W^|U=o W < C w \\H\\ w i,oo ix) , (12) 

l<i,j<m 

\\W-WWoo~mzx \\W ij -W! j \\ L ~(X)<C! / r\\H-H'\\ w i,- m , (13) 

whenever two data sets H and H' jointly satisfy condition (7) with the same triple 
(C,r, Co). After proving solvability of this system, one is able to reconstruct the 
functions Ui and then to reconstruct a as described below. Uniqueness and stability 
of the solution to (11) with respect to the drift fields Wij relies on the fact that 
— Cp 1 = — ((n — 2)a + l) is not an eigenvalue of the operator Pyy '■ W — > H defined 

by 

P w ■ v i-> P ff v = [P w w] iei = A fl ] (% • V Vj ) e u (14) 

where A^, 1 denotes the inverse of the Dirichlet Laplacian on X , and where we have 
defined the space H := {Hl(X)] m , which makes (H, || • ||-«) Hilbert once equipped 
with the norm 

m m » 

= X>i|&i =J2 / l V ^! 2 dX > V =(«!,..., « m ). (15) 
i=l (=1 Jx 

When the coefficients Wij are bounded, we show that the operator Pw is compact 
and its operator norm satisfies the estimate ||Pw|l — V™W II Halloo (see lemma 
5.1), where HA^ 1 )! denotes the operator norm of A^, 1 : L 2 {X) — > Hq(X). As 
a consequence, the system (11) satisfies a Fredholm alternative which will provide 
uniqueness and stability as stated in the next proposition, for all a € M when n = 2, 
and for all a but possibly a discrete set (possibly converging to — (n — 2)" 1 ) in the 



interval 



whenever n > 3. 



Proposition 1 (Stability of the strongly coupled elliptic system). Let vector fields 
{Wij, W|j}i<jj<ro belong to L°°(A) and such that — Cp 1 is an eigenvalue of neither 
Pw nor Pft"- Let u, u' be the unique solutions to (11) with same illumination g 
and respective drift terms W , W' . Then we have that u — u' S T-L and satisfies the 
stability estimate 

||u-u'|| w <C||W-W'||oo. (16) 

Remark 3. In the case n = 2 or (a — with m = n), we can recast (11) as a 
coercive system in divergence form, the injectivity of which follows immediately. 
These cases correspond to cf = 1- 

Once the solutions u, arc reconstructed, one may reconstruct a using a formula of 
the form a = -ffii/|Vwi| 2 . However, such a formula may not offer the best stability 
estimates. Another reconstruction strategy is deduced from (9), which can be recast 
locally as 

V(<r- 2 «) = -^f;(V(Dir ij ).V« i )V« i , a^O, (n-2)a^-l, 

(17) 

Vloga = - J2 (V • VuiJVu,-, a = 0, 
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where Vui, . . . , Vw n denote the n linearly independent gradients. As in the ODE- 
based reconstruction procedure, we can devise an ODE-based algorithm to recon- 
struct a locally from formula (17). We then arrive at the following stability result. 

Theorem 2.4. Let the conditions of Proposition 1 be satisfied. Then the corre- 
sponding a, a' satisfy the estimate 

\W~ 2a - a'- 2a \\ m(x) <C\\H-H'\\ w ,,oa {x) , a^O, (n - 2)a ? -1, 
|| logcr - logcr'|| ff i (x) < C\\H - H'\\ w i,^ {x) , a = 0. 

Note that a necessary and sufficient condition for the unique solvability of (17) 
on a simply connected domain is that the exterior derivative of the right-hand 
side (seen as a one-form) vanish by an application of the Poincare lemma. The 
compatibility condition that arises here takes the form of a quadratic equation in 
the components of Vu,; that is difficult to ensure as it depends on the unknowns. 



3. Geometric setting and proofs of lemmas 2.2 and 2.1. Defining geometric 
notation for now, let us first denote the Euclidean orthonormal frame e.; = d x % 
and e l = dx % . For < k < n, A k denotes the space of k— forms. We recall the 
definition of the Hodge star operator * : A fe — s- A"~ fe for < k < n, such that for 
any elementary fc-form dx 1 = dx 11 A • • • A dx %k with I = . . . , if-), we have 

-kdx 1 — adx J , where a = sign (1 . . .n) <— ¥ (/, J). (19) 

We recall the following useful identities, see e.g., [23]: 

**=(_!)*(«-*) on A fc , *(u b A*v b ) = u-v, *d*u b =V-w, u.wGA 1 . 

We now prove Lemma 2.2 which is the cornerstone of our explicit reconstructions. 

Proof of Lemma 2. 2. Because S\ (x) , . . . , S n (%) is a basis of K ra at any point i£X, 
a vector V can be represented in this basis by the following representation (x is 
implicit here) 

V = H ij (V ■ Si)Sj. (20) 

For j = 1, . . . , n, let us introduce the following 1-forms: 

X) := (-l) n -\ * (Si A-ASfj, (h,...,i n . 1 ) = (l,...,j,...,n), (21) 

where the hat indicates an omission and crj = (— l)- 7 " 1 is the signature of the 
permutation (1,2 ... ,n) i— > (j, 1, . . . , j — l,j + 1, . . ., n). At each x € fl, the vector 
Xj(x) obtained from Xj(x) by "raising an index" can also be seen as the unique 
vector obtained by the Riesz representation lemma that corresponds to the linear 
form Dj : E" -> E such that for any V £ R n , 



Dj(V) = det(Si(i), . . . , Sj_i(aO, V, S j+1 (x), ...S n (x))= X 3 (x) ■ V. 
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We now show that the vector fields Xj satisfy a simple divergence equation. We 
compute 

V • X S = *tf * X) = a, * d(Sl A • • • A Sl_ t ) 

n-1 

= a, A ' ' ' A d5 t A ■ ■ ■ A 

fe=i 

n-1 

= a, * 2 (-l) fe ^ A • • • A a(F b A £ ■ J A • • • A 

fe=i 

= {n- l)a*(F b A*X b ), 

and using the identity *(w b A *v b ) = u ■ v, we deduce 

V ■ X,- = (n - l)aF ■ X,-, j = 1 . . . n. (22) 

The decomposition of Xj in the basis S\ , . . . , S n may be obtained by computing its 
dotproducts with S±, . . . , S n . Indeed, for k j, there is an I such that i\ = k and 
we have 

Xj • Sk = det(5i, . . . , Sj-i, Sk, Sj+i, . . . , S n ) = 0, 
by repetition of the term Sk in the determinant. Now if k = j, we have 

Xj ■ Sj = det(Si, ...,S n )=detS = D. 
Using formula (20), we deduce that Xj admits the expression 

A, = l)U ,J S,. 

Plugging this expression into equation (22), and using V • (ipV) = Vip ■ V + ipV ■ V, 
we obtain 

V{DH lJ ) ■ S t + DH ij \> ■ Si = (n - l)aF ■ {DH V] S l ) 
<=> V(DH ij ) ■ Si - DH ij (l - a)F ■ S d = (n - l)aDH ij F ■ S, 
& V(DH ij ) ■ Si = r,.. '/;//"/• • Si. 
Finally using the representation (20) for F itself yields 

F = (H ij F ■ Si)Sj = ^(V(DH ij ) ■ Si)Sj. (23) 
We can also recast the previous expression as follows 

F = c F [H ij {VlogD ■ S t )S 3 + ((ViT 3 ') • St)Sj] = c F [VlogD + {{VH ij ) ■ S^Sj] , 

(24) 

and the proof is complete. □ 

We now give a proof of Lemma 2.1, which guarantees the existence of illumi- 
nations that ensure condition (7) and thus justifies the two global reconstruction 
approaches. The CGO constructions, introduced in [7] in this context, generalize 
those defined in [5]. 

Proof of lemma 2. 1 . Since a is bounded from above and below by positive con- 
stants, it suffices to study the case a = \ since we have for any a\,a^ € K, 

dct(cr ai Vlii, . . . , (T ai Vll n ) = o* (ai_aa) dct(cr Q2 V«l, • ■ • , CT" 2 Vli n ). 
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Consider the problem V • <j(x)Vu = on W 1 with a(x) extended in a continuous 
manner outside of X and such that a equals 1 outside of a large ball. The con- 
struction requires sufficient smoothness of a in order to be valid. Let q(x) = — A ^ 
on R". We assume that q G H^ +1+e (M n ), which holds if a - 1 G iff + 3+£ (R n ) for 
some e > 0, i.e., the original <7|x G H% +3+£ (X). Note that by Sobolev imbedding, 
a is of class C 3 (X) while g is of class C 1 (X). With the above hypotheses, we can 
apply [7, Corollary 3.2] which states the following. 

Let v = y/au so that Av+qv = on R n . Let p G C™ be of the form p = ^(k+ik- 1 -) 
with k,k x G S n_1 , k • k- 1 = 0, and p = \p\/y/2 > 0. Thus, p satisfies p ■ p = 
and e p :c is a harmonic complex plane wave (hence the name of complex geometrical 
optics solutions). Now, it is shown in [7], following works in [8, 22], that 

v p = V^u p = e" x (l + ib p ), pib Plx = 0(1) in C l (%), 

with (A + q)Vp = and hence V • <tVu p = in W 1 . We have used again the Sobolev 
imbedding stating that functions in H% +k+£ (Y) are of class C k (Y) for a bounded 
domain Y. Taking gradients of the previous equation and rearranging terms, we 
obtain that 

yfaVup = e p ' x (p+(p p ), with ip p := V^ P + %b p p - (1 + ip p )Vy/a. 

Because V-y/cr is bounded and pip P i X = m ^POi the C n -valued function 

ip p satisfies sup^|<p p | < C independent of p. Moreover, the constant C is in fact 
independent of a provided that the norm of the latter is bounded by a uniform 
constant in H% +3+£ (X). 

Both the real and imaginary parts of u p , denoted up* and , count as solutions 
of the free-space conductivity equation, thus \/aVu^ and -y/erVitp can serve as 
vectors Si. More precisely, we have 

v^V< - pe*" ((k + p- 1 ^) cos(pk^ • a:) - (k 1 + p~ l ^ p ) sinfck 1 • x)) , 
VSVuj* = pe pk ' x ({k 1 - + p- Vp) cos(pk x ■ x) + (k + p~ V* ) sin(pk x • x)) . 
Case n even . Set n = 2p, define p x = p(eii + for 1 < Z < p, and construct 

S^i-i = VffVw* and S* 2/ = VaVu* , 1 < Z < p. 
Using (25), we obtain that 

det(Si, . . . , S n ) = p n e 2 "^ ** (1 + f(x)), 

where linip^oo sup^ |/| = 0. Letting p so large that sup^ |/| < | and denoting 
7o := mm xeX (p n e 2 P^=i X2 ') > 0, we have inf^ det(Si, . . . , S n ) > ^ > 0. We 
conclude after the next paragraph. 

Case n odd . Set n = 2p — 1, define p ; = p(e2; + ie2(-i) for 1 < I < p — 1, and 
p p = p(e n + iei) and construct 

S-w-i = VoVup, and 5*2/ = v'oVup, - 1 < * < J>- 
Using (25), we obtain that 

det(5i, . . . , S n -!, S n ) = p n e"^+ 2 ^ x ^ (- cos(pzi) + h(x)) , 
det(Si, . . . , S n -i, S n+1 ) = p n e"^+ 2 ^ x ^ (- sm(p Xl ) + f 2 (x)) , 
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where lim p _>oo supy |/,-| = for i = 1, 2. Letting p so large that supy(|/i|, I/2I) < \ 
and denoting 71 :~ Ttun xt --^(p n e p ^ Xn+2 Sf=i > q. we have that 



|det(Si,...,S n _i,5 B )| > .TeA'n m e y,^ 



J717T 



|det(5i,...,5„_i,5' n+ i)| > -y-, ieln Ipxi e [ ) +m7r 



4 ' t \6 6 

where m is a signed integer. Since the previous sets are open and a finite number of 
them covers X (because X is bounded and p is finite) , we therefore have fulfilled the 
desired requirements of the construction. Upon changing the sign of S n or S n +i on 
each of these sets if necessary, we can assume that the determinants are all positive. 
Conclusion . In each of the previous cases, let {<7i}i<;< m be the traces of the solutions 
defined above with m = 2[_ Ii ^-J . These illuminations generate solutions that satisfy 
the desired properties of maximal rank and positive determinants. By continuity 
arguments, any boundary conditions gi in an open set sufficiently close to gi will 
ensure that the maximum of the determinants stay bounded from below by Co > 0. 
This concludes the proof of the lemma. □ 

4. The ODE-based method. In this section, we extend the results presented in 
[5] to general dimension and for a more general class of measurements (described by 
the coefficient a). We first need to introduce standard geometric notation, without 
which the derivations become quickly intractable. 

4.1. Definitions, notation and identities. We work on a convex set fi C K n 
with the Euclidean metric g(X,Y) = X Y = SijX'Y^ on R n . Following [15], we 
denote by V the Euclidean connection, i.e. the unique connection that is torsion- 
free, and compatible with the Euclidean metric in the sense that 

V X (Y ■ Z) = (V X Y) -Z + Y- (V X Z), 

for smooth vector fields X, Y, Z. On zero- and one-forms, this connection takes the 
expression: 

V x f = X-Vf = X i d i f, and V X Y = (X ■ VY^ej = X^Y^ej, 

for given vector fields X = X % e,i and Y = Y z ei. An important identity for the 
sequel is the following characterization of the exterior derivative of a one-form to 

du(X,Y)=W x (u(Y))-V Y (w(X))-w([X,Y]), (26) 

or equivalently in the Euclidean metric, writing lo = Z~° for some vector field 

Z-[X,Y] = V x {Z-Y)-V Y {Z-X)-dZ\X,Y), (27) 

where the Lie bracket (commutator) of X and Y coincides with (and thus may be 
"defined" here as) [X, Y] = V xY — VyX by virtue of the torsion-free property. 

A frame refers to an oriented family E = (Ei, . . . , E n ) of n vector fields over Q 
such that for every x € 0, (Ei(x), . . . E n (x)) is a basis of T x Vl = K™. For a given 
frame E, we define the Christoffel symbols (of the second kind) with respect to this 
frame, by the relations 

Ve.Ej = E k , i.e. r?, =g^V Ei E } ■ ■ E p , where 

(28) 

g ij = E i -Ej and g pq = (g ) pq . 
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The following very useful identity allows us to compute the Christoffcl symbols from 
inner products and Lie brackets of a given frame (see e.g. [15, Eq. 5.1 p. 69]): 

2(V X Y) ■ Z = V X (Y ■ Z) + V Y {Z ■ X) - V Z (X ■ Y) 

- Y ■ [X, Z] - Z ■ [Y, X] + X ■ [Z, Y], 

where X, Y, Z are smooth vector fields. 

For a vector X = X J ej, we want to form the matrix of partial derivatives 
(djX l )ij. Geometrically, gradients generalize to tensors via the total covariant 
derivative, which maps a vector field X to a tensor of type (1,1) defined by 

VI(w,y)=w(VyI). (30) 

In a given frame E, we may express VEi in the basis {Ej ®^}^ fc=1 of such tensors 
by writing V-Ej = aijkEj ® E k and identifying the coefficients aij k by writing 

VEi(E b p ,E q ) = E b p (V E q Ei) = V E q E, ■ E p = g p7 T r qi , 

and also 

V Ei(Ep, E q ) = a^kEj ® E b k (E p ,E q ) = aij k gj P g kq . 

Equating the two, we obtain the representation 

V-E, = g^Ej ®E b k = ff «V p (V B ,^ • E P )E, ® E{. (31) 

The theory of the following sections proves that all partial derivatives of a frame 
(given in (31)) are uniquely determined by inner products gij and by Lie brackets, 
as (29) indicates, or equivalently by exterior derivatives, as (27) expresses. These 
derivations will be carried out first for the S frame and second for the R frame with 
values in the space of rotations SO(n, ffi). 

4.2. The S frame. We now study the properties of the S frame. S is a frame 
provided that the determinant condition inf^go det S > cq > holds. Our objective 
in this section is to find an expression for VSi that allows us to solve for Si by the 
method of characteristics. We have seen in the preceding section that this involved 
calculating the Lie brackets (commutators) of the vectors composing the frame. For 
I < i < j < n, we have 

[S i ,S j ] = H kl ([S i ,S j }-S k )S l . (32) 

Now using (27) we write 

S k ■ [S„ Sj] = V Si (S k ■ Sj) - V Sj (S k ■ S,) - dSliS,, Sj) 

= S, ■ VH kj - Sj ■ VH ki - aF b A S^Sj) 

= Si ■ VII,, - Sj ■ V//,,, + n,: //;. ,/•• • Si + H M F ■ Sj). 

Plugging this into (32) and using that H kl Hkj = %, we obtain the Lie brackets 
[Si,Sj] for 1 < i < j < n : 

[S u S 3 ] = H kl [WH jk ■ Si - WH ik ■ SjjSt + a((F ■ S j )S i - (F ■ SJSj). (33) 

Returning to the computation of VSi using (31), we combine (33) with (29) to 
arrive at 

2(Vs„£i) • S p = V.s q Hip + VstHpq — V s p Hqi 

Si • [Sq : Sp] Sp • [Si : Sq] -\- Sq • [Sp, Si\ 

= VH iq ■ S p + VH ip ■ S q - VHpq ■ S t + 2a{H pq {F ■ Si) - H qi (F ■ S p )). 
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Plugging this expression into (31) (expressed in the S frame), and using H^W k = 
Sik, we obtain 

2WSi = 2m k H^{v Sq s t ■ s p )Sj ® si 

= //''' // '•"( V//,, • S p + VH ip ■ S q - VHpg ■ Si 
+ 2a(H pq (F ■ Si) - H qi {F ■ S P )))S 3 <g> S b k 

= (H jp u ik ■ s p + m l r n ■ s q + v// 4 • s, 

+ 2a(W k (F ■ Si) - H™S ik (F ■ S P )))S 3 ® Si 

where we have used ViP = —Hi p (V H pq )H qk and have defined 

U jk := (VH jp )H pk = -H Jp VHP k , l<j,k< n. (34) 

Using formulas W k S 3 ®S{ = l n := e i (g)e i and H kl (V-S k )S l = V for any smooth 
vector field V, we obtain for 1 < i < n 

VSi = i (U ik ®Si + S k ®U ti lk + (VF t: • Si)Sj ® Sty + a{F ■ Si)!* - aF ® S\. 

(35) 

Using (23), we observe that VS; is equal to a polynomial of degree at most three 
in the frame 5* with coefficients involving the known inner products Hij . For each 
1 < i,k < n, d k Si is nothing but V e( .S; = VSi(-,efe), which can be obtained 
from (35). Denoting S := [Si \ . . . \ S n ], we are then able to construct the system of 
equations 

2 

n 

d k S=J2Q k p S ^ S"=JJsf, l<k<n, (36) 

l/3|<3 i=l 

where Q k p depends only on the data and is an n 2 -indcx. This redundant system 
can then be integrated along any curve (where it becomes a system of ordinary 
differential equations with Lipschitz right-hand sides ensuring uniqueness of the 
solution) in order to solve for the matrix- valued function S. 

4.3. The orthonormal R frame. The above system (36) involves a priori n 2 
unknowns since the matrix S does not necessarily have any useful symmetries. 
However, we know the inner products H = S T S, i.e., a matrix of dimension \n{n + 
1). We therefore hope to be able to find a closed-form system involving |n(n — 1) 
dimensions. This is the dimension of the orthonormal R frame. 

We now provide the details of remark 1. From the frame S, we build an oriented 
orthonormal frame R = [Ri \ ■ ■ ■ \R n ] (or equivalently, an SO n (R)-valued function) 
from a matrix- valued function T(x) = {Uj{x)}i<i,j<n that satisfies the relations 
T T T = H^ 1 and detT > at every x € f2, as well as a stability property of the 
form 

||T - T'\\ w ^ {x) <C T \\H- H'\\ w ^ {x) , (37) 

where Ct > depends only on the way we construct T from H . T can either be 
constructed by the GS procedure or by setting T = H~?, the positive square root 
of H . The stability statement (37), first proved in the GS case for n — 2, 3 in [5], 
can be obtained for both GS and T = H~?, see [16] for proofs of these statements. 
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The function R := ST T satisfies everywhere R T R = I„ and deti? = f, hence R 
is an SO n (R)-valued function. The column vectors of S and R transform according 
to: 

Rt = Uj Sj , S, /"/.*,. i = 1 . . . n. (38) 
We also define for 1 < i, k < n 

V lk := (VUj)t jk , Vf k :=±{V ik + V ki ) and V? k := ^(V ik - V ki ). (39) 

We are brief on the derivation of the gradient system for R as it is very similar to that 
of the S frame. The system of equations (4)- (5) together with the transformation 
rules (38) allow us to derive the following system of equations for the R frame: 

V • Ri = V lk ■ R k - (1 - a)F ■ R u (40) 

dR\ = V> k A R k + aF b A R\, 1 < i < n. (41) 

From this system, we express F in the R frame as 

F = c F (V logD + ((Vij + Va) ■ R i )R j ) • (42) 

Equation (42) can also be derived directly from (9) and the transformation rules 
(38). Then, using equation (41) and formula (27), the Lie brackets [Ri, , Rj] of the 
vectors take the form, for 1 < i < j < n: 

[Ri,Rj] = (-Vpj ■ Ri + V pi ■ Rj)Rp + a((F ■ Rj)Ri - (F ■ R l )R 3 ). (43) 

From (43) we deduce the Christoffcl symbols relative to the R frame: 

r|. = Vf k ■ Ri + V? k ■ Rj - Vij ■ R k + a(F ■ R 3 )S lk - a(F ■ R k )8 l0 . (44) 

Finally, in the orthonormal case, the expression of the gradient reduces to V 'Ri = 
T ki Rj ® R b k , from which we deduce that 

Vi? 4 = R k ® - V& ®Ri + (Vf k ■ Ri)Rj ®R b k + a(F ■ Ri)I n - aF ® R\. (45) 

As for the S frame, the R.H.S. of (45) depends polynomially on R and on the data. 
This system can thus be solved for the vectors Ri via ODE integration along any 
curve in a connected domain and provided that we know the R frame at one point. 
In practice, this system is less expensive to integrate than (36) since the R frame 
can be locally parameterized with n(n — l)/2 scalar functions (such as the Euler 
angles) whereas the S frame requires n 2 scalar functions. 

4.4. Global reconstruction algorithm. The proof of the stability theorem 2.3 
can be found in [5] in dimension n = 3 with a — | (although the proof would be 
identical in arbitrary dimension). In that paper, the theorem is proved using the 
system for the rotation matrix R and thus requires the extra stability condition 
(37). This condition is necessary only if we reconstruct a via the R frame. The 
same stability result can be obtained without this requirement if we reconstruct 
a via the S frame directly. In the latter setting, the proof is quite similar to the 
one in [5] with the further simplification that we do not need to change bases when 
switching sub-domain fi^. The system of ODEs that one must solve based on the 
gradient system (35) is well-posed since the function S satisfies a priori the uniform 
bound 

m 

\S(x)\ 2 =Y, H » <™\\H\\co, 

i=l 
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and the right-hand side of (36) is Lipschitz in S over the set {S : X ->• M. nm , ||S||oo < 
\J m\\H \\oo} as a polynomial of the components of S, and using the fact that the 
polynomial Q k are bounded; see [5] for additional details, which we do not reproduce 
here. 

5. The elliptic method. 

5.1. Derivation of system (11) and equivalent formulations. 

5.1.1. The case m = n. In this case, condition (7) is satisfied with the partition 
O = {X}, N = 1. Equation (23) can be rewritten as 

VlogCT = j^<j 2a (V(DH kl )-Vu k )Vui. (46) 

Rewriting the conductivity equation (1) as 

Aui + V log a ■ Vu-i = 0, 
and plugging (46) into it yields the coupled elliptic system of equations 

= Au t + ^(V{DH kl ) ■ Vu k )a 2a Vui ■ Vi k = Au, + c F W lk ■ Vu k , (47) 
where we have have defined 

Wtk := ^-V{DH lk ) =VlogDS ik + H a VH kl , 1 < i,k < n. (48) 

From the last form of Wi k , we derive (12) and (13) since the denominators only 
involve D which is bounded away from zero and the rest is polynomial in the Hij 's 
and their derivatives. 

Multiplying (47) by DH pt and writing it in divergence form, one obtain the 
following equivalent formulation to (47) in variational form: 

-V • {DH^Vui) + (1 - c F )V{DH pl ) ■ Vui = 0, 1 < p < n. (49) 

5.1.2. The case m > n. In the case where we have m > n solutions, we can still 
define V log a over the entire domain X using a partition of unity that is subordinate 
to the open cover 0, call it {Lpi}fLi- Then we can define V log a globally over X by 
writing 

N 

Vlogcr = Vlogp-ln^, 

where the restrictions are constructed from the n solutions of positive determinant 
on each f2^. Each of these restrictions can still be written in the form 

m 

Vlogoisi, = c F a 2a ^2 ( F 3k\ui ■ V%)Vitfe, where 
j,k=i 

{0 if j £ r(i) or fc ^ r(i) 

(irWll^- 1 ) if (j, k) = (r(i) OJ r(i) b ), 

with D T{ ^ = Vdct W®. Thus we can patch these formulas together into a globally 
defined 

m N 

V log o := c F a 2a (F jk ■ Vuj)Vu k , where F jk = ^ F jk \ Sh ip t . 

j,k=l i=l 
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Plugging this expression into the conductivity equation yields the coupled elliptic 
system 

= Aui + V log a ■ Viti = Aui + cwHikF^v. ■ Vu,-. 

(50) 

Ui\dx = gi, l < t < m, 

and one arrives at a system of the form (11) by setting Wij := H^Fj^ for every 
1 < i,j < m - In this case, the stability inequalities (12) and (13) can be derived 
using the fact that 

||VF|| L oc (x) < i muc r ||W||L-(n < ), 

and noticing that on each Oj, Wij is either zero or locally defined by (48) (and 
weighed by tpi) whose expression has been proved to be stable. A similar argument 
holds for proving (13) thanks to the fact that the partition of unity {^Pi} is the same 
for two data sets H,H' that jointly satisfy (7) with the same triple (O, r, c$). 

5.2. Uniqueness and stability results. 

5.2.1. Proofs of Proposition 1 and Theorem 2.4- Let us assume a system of the form 
(50), where the vector fields Wij belong to L°°(X). Assuming the illumination g to 
be in [Hi (dX)] m , we use a lifting operator to define functions {u^}™ 1 E [H 1 (X)] m 
of traces g at dX. Defining the unknown Vi = Ui — Wi, we are now left with analyzing 
the solvability of the system 

Av t + c F Wij ■ Vvj = hi (X), Vi\ dx =0, 1 < i < m, (51) 

where hi := Aw t + c F Wij ■ Vwj e H^ 1 (X), (52) 

as well as the stability of its solution with respect to the vector fields Wij. H^ 1 (X) 
denotes the dual space of Hq (X). 

As described in section 2, we apply the inverse of the Dirichlet Laplacian to (51) 
and obtain the system of integral equations 

Vi + c F A' D 1 {W l j ■ Vvj) = A^hi, 1 < i < to, 

which can be recast in vector notation as 
(I + cj?P»-)v = f , where 

_i _i (53) 

P w v := [P w v]iej = A D (Wij ■ Vvjfe, and f := A D hi e,. 

Because A^ 1 is continuous in the functional setting H~ 1 (X) — > Hq(X) (see [11]), 
it is clear that f belongs to the space T-L. We now have the following: 

Lemma 5.1. Assuming that the vector fields Wij G L°°(X), the operator Pyy : 
% — > % defined in (53) is compact, and its norm satisfies 

\\P W \\ < vMIA^ 1 1| H^Hoo, ||W||oo = max H^IU- (54) 

l<i,j<?n 

Proof. As can be seen in [11] for instance, the operator A^ 1 : L 2 (X) —> H 2 (X) 
is bounded. Therefore, by the Rellich compactness theorem, the operator A^, 1 : 
L 2 (X) — > Hq(X) is compact and of norm denoted by ||A^ ||. Now P is also 
compact since each of its components is the composition of the continuous operator 
W^ ■ Vvj e L 2 (X) with the compact operator A D l : L 2 (X) -> H^(X). 
Moreover, for v£?i and every 1 < i,j < m, we have the obvious bounds 

\\A^(W l3 ■ V Vj )\\ H i < IIA^IHI^ • V Vj \\ L 2 < IIA^IIII^IUII^II^, 



16 



FRANQOIS MONARD AND GUILLAUME BAL 



and thus 

\\[Pwv]i\\% h <\\A„Y\\W\\Uv\\ 2 n. 
Summing over i proves (54). The proof is complete. □ 

As a consequence of lemma 5.1 and by virtue of standard compact operator 
theory (e.g. [11, Theorem 6 p 643]), we have the following facts: 

• is eigenvalue of Pw, which corresponds to the case (n — 2)a = — 1, a value 
for a that we exclude from our analysis, 

• the remaining spectrum of Pw is point spectrum and consists of at most a 
discrete sequence of values that is either finite or converges to zero. 

Finally, the operator I + c F Pw G satisfies a Fredholm alternative. Therefore 

it suffices that — Cp 1 ^ sp (Pw) in order to obtain uniqueness and stability of 
the solution of (53) and therefore of the solution of (11) as well. The proof of 
Proposition 1 makes these statements more precise. 

Proof of Proposition 1. Let W, W have their coefficients in L°°(X) and such that 
— Cp 1 ^ sp (Pw)Usp (P w ,), an d let v, v' e H solve the system (51) with respective 
drift terms W, W and same illumination g. By virtue of the Fredholm alternative, 
the operators I + c F P\y and I + c F Pw' are invertible with continuous inverses in 
C(H). 

Applying the inverse Dirichlet Laplacian to both systems, we obtain the systems 
(I + c F P w )v = i, and (I + cpPwW = f'- 
Taking the difference of both systems, the resulting system reads 

(I + c F P w )(v - v') = f - f - c F Pw-w'v'. (55) 
The first difference in the right-hand side of (55) may be bounded by 

||f - f'\\n = Wcf^d 1 {{Wn -Wlj)- Vtwj) ei\\ n < Cc F \\W - VF'IU, 

where the constant C depends on ||A^ [l^^-i #i) and maxKj< m llSill^i^^- ^-P" 
plying lemma 5.1 to the operator Pw-w- the second difference in the right-hand 
side of (55) may be bounded by 

\\P W -w>v'\\n < v^llW-^'lUIIA^IHIv'llw 

< V^\\W- W'HoollA^IIIKH-CFPwO^llllflk- 

Combining the last two estimates with (55) and the fact that I + c F Pw is invertible 
with continuous inverse in C{H), we arrive at 

l|v- v'Wu < c"||(i + c^Pw)- 1 !!!!^ - v^'Hoo, 

for some constant C > 0. Since u u' = v v', this concludes the proof. □ 

We now conclude with the proof of theorem 2.4. 

Proof of theorem 2.4- We focus on the case a/0 and (n — 2)a ^ —1. The proof 
for a = is identical up to small changes in notation. Let H, H' have their compo- 
nents in W 1,oa (X) and jointly satisfy (7) with the same triple (0,t, cq). Then the 
families of vector fields W and W' have their coefficients in L°°(X) and we further 
assume that — Cp 1 £ sp (Pv^)Usp (P w ,). Let v,v' e % solve the system (51) with 
respective drift terms W and W and same illumination g, and let a, a' be the cor- 
responding conductivities. Without loss of generality, we work on one of the open 
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sets Qi £ O and renumber the n solutions whose gradients are linearly independent 
from 1 to n. The result will then hold provided that we have X C Uji^fij and thus 



N 

I -2a r-2a\\2 ^\^II„-2q i-2a\\2 



For a given Hi £ O, and defining Vij := —2a £ £-W(DH' l: '), we write, using equality 
(17) 

V( ( t- 2q - a'- 2a ) = ((Vy - ^) • VuOV^ + (F/. • V( Ul - <))V« 3 - 

+ (^-v«J)vk--u^). 

Similarly to the vector fields Wij (48), the vector fields Vij satisfy estimates of the 
form 

llV^-lloo < CVH-ffllwri— and WVij-V^WooK&yWH-H'Wwi.*, l<i,j<n. 

(57) 

Since H is bounded and a, a' are assumed to be bounded from below by a constant 
Co > 0, each of the Vit^ is uniformly bounded by y Ha / a Q < *J\\H Hoo/^o- Taking 
L 2 norms over and using the triangle inequality in (56), we obtain 

HV^- 2 " - a'- 2a )\\ L2{Ui) < IIVttilUHVuillooll^ - Vl 3 \\ L , 

+ 2||^.|| 00 ||V(« i -t4)|| ia ||V« j || OOJ 

which by virtue of proposition 1 and estimates (57) yields an estimate of the form 

||V(a- 2Q - a'- 2Q )|| L2(ni) <C\\H- H'\\ w ^ {x) . (58) 

Further, from the pointwise relations a~ 2a = |V«i| 2 / H\\ and similarly for a', we 
write 

-2a '-2a _ * V7f„, „/ \ V7f„, i \ i (V IT' \ \^ U ll\ 2 



a — a 



-v(«i - o ■ v(«i + o + (h u - frd: 



-"11 ilnii 11 

Taking L 2 norms and using the triangle inequality, we obtain that 

\W~ 2a - c7'- 2a \\ L 2 {ni) < (||v(«i - uDMiivuiiioo + livuiiu) 

+ l|H 11 -H 11 y 00 ||F 1 - 1 1 l| 00 ||Vt i / 1 |U 2 ||Vu / 1 || 00 )||Hr 1 1 || 00! 

which again yields an estimate of the form 

\W~ 2a - «'- 2a hHn,) <C\\H- H'\\ w ^ [x) . (59) 
Combining (58) and (59), we arrive at 

h~ 2a - a f - 2a \\m(fl t ) <C\\H- H'\\ w ^ {x) , 
for every Hi £ O. This concludes the proof. □ 



5.2.2. Discussion on the spectrum ofP\y- 
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The two-dimensional case . In the case n = 2, [2. Theorem 4] guarantees that we 
can pick m = n. In so doing and using the form (49) of the elliptic system together 
with the fact that cp = 1 for all a £ R, we arrive at the system 

V • (DH^Vm) = 0, u p \ dX =g p , p = l,2. 

The weak formulation of the corresponding problem with homogeneous Dirichlct 
conditions involves the bilinear form 

B(v,v) = / DH pt \7 Vl ■ \7v p dx. 

Since if -1 is uniformly elliptic over X and infx D > co, this bilinear form is coercive 
over T-L as seen from the following calculation 




2 



where Am stands for the largest eigenvalue of H, for which we have, pointwise (xm 
designates a unit eigenvector associated with Am) 

Am = {xmiHx m ) = HijX M>i XM,j < \\H\\oJj^2x 2 Ml + x 2 M j = nWH^, 

and hence the estimate 

B(v,v) = / DH*Vvi ■ Vv p dx > CD(n||fr|| 00 )- 1 ||v|&. 
Jn 

Therefore by virtue of the Lax-Milgram theorem, the system (47) admits a unique 
solution in %. In particular, this shows that — Cp 1 is not an eigenvalue of Pvy in 
this case for any a £ R. 

The case n > 3 . Using the fact that the spectrum of ||PvHI i s bounded in norm 
by ||Pjj-|| and that Pu' is compact, the elliptic system admits a unique and sta- 
ble solution, except for a discrete set of values — Cp 1 £ [— ||PvHU|PvHI] possibly 
converging to zero. In terms of a, this corresponds to almost all values of a £ R ex- 
cept a sequence {ak} taking values in the interval [ ^JjL"" 1 ] and possibly 
converging to — (n — 2) _1 . 

The special case a — 0. In this case we have cf = 1. This implies that whenever 
one can ensure the positivity condition (7) with only m = n solutions (e.g. in even 
dimension and using lemma 2.1), one can rewrite the system into the form (49) with 
term 1 — cp = 0, that is 

V • {DH pi Vui) = 0, 1 < p < n. 

Using the same arguments as in the two-dimensional case, this system is coercive 

and therefore ensures that —1 is not an eigenvalue of P^. 

Conclusion . As a conclusion of this discussion, the following statements hold: 

1. if n = 2, we have sp P w n {-((n - 2) a + 1),q£I] = 0, 

2. if n > 3, then sp Pw H{— ((n — 2)a+ 1), a £ M.} consists of at most a sequence 
{-((n - 2)a k + 1), k = 1, 2 . . .} where a k belongs to [ ""^J" 1 , l|P ^~ 1 ] and 
possibly converges to — (n — 2) _1 . In the case m = n, the value is excluded 
from the latter interval. 
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6. Constraints, reconstructions, and compatibility conditions. The ODE- 
based reconstructions use the full redundancy of the data to construct an overde- 
termined system of equations for the vectors Si (or Ri) and the vector V logo - . The 
PDE-bascd method defines a well-posed system of equations for the scalar quanti- 
ties Ui and an overdetermined system of equations for vector V log a. Each of these 
ovcrdctermined systems needs to satisfy compatibility conditions in order to admit 
a solution. In this section, we aim to extract information from the over-determinacy 
of the system. We first revisit the two-dimensional case and use the redundancy to 
extract explicit reconstruction algorithms in the setting a / ^. We then consider 
the case of arbitrary dimension and show that the compatibility conditions that 
data must satisfy in order for the aforementioned systems to have solutions take 
the form of vanishing appropriately defined curvatures together with the cancella- 
tion of a given two-form. These conditions generate quadratic functionals of the 
unknown vectors Si or Ri whose pratical applicability is discussed below. 

6.1. Reconstructions in two dimensions. In this section, we revisit the two- 
dimensional case which was first solved in [5, 10] and generalize the approach to 
the case a / ^. In that approach, the reconstruction of F = V log a requires the 
reconstruction of a function 9 : X — > S 1 that characterizes the unknown information 
about the frames S or R. We consider the S^W-valued R = (Ri, R2) frame and 
parameterize it as R\ = (cos8,s'm8) T and R2 = JRi, with J := [ ° With the 
notation := Ri ® Rj for i,j = 1, 2, we recast equations (42) and (40) for a e R 
as follows: 

2 

V log a = F = V log D + 2 ®ij V ij > ( 60 ) 

*ij'=i 

V ■ Ri ; = V lk ■ R k - (1 - a)F -Ri, 1 = 1, 2. (61) 

We next derive an equation for V#, which by construction is nothing but [R 2 , Ri]. 
We have: 

V0 = [R 2 ,Ri] = v R2 Ri - v Rl i? 2 = -rl 2 R 2 + r^fli, (62) 

where the Christoffel symbols 1 and T 22 are given by 
T 2 n = V Rl R! ■ R 2 = [V, = -V • R 2 

= -((2a - 1)V 22 - (1 - a)N) ■ R 2 - ((2a - l)V& - V£) ■ Ri, ^ 

. (63) 

T\ 2 = V fl2 i? 2 • R x = -[V,R 2 ] = -V • Rt 

= -((2a - 1)V U - (1 - a)N) ■ Ri - ((2a - l)V? 2 + V? 2 ) ■ R 2 . 
By orthonormality the other Christoffel symbols are given by 

r n = r 21 = V\ 2 = T 22 = 0, r} 2 = — r lx , and r 2 i = —re- 
plugging the expressions (63) into (62), we arrive at 

V8 = V& - (1 - a) JVlogD + (2a - l)($ 2 i^ii - $12^22 + ($22 - $11)^12)- 

(64) 

Using the following identity 

$ 21 Vn - $12^22 + ($22 - $11)^2 = J($nV n + $ 22 V 22 + ($12 + $21)^2) 

= ij( j p-Vlog£»), 
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equation (64) may be recast as 

V9 = V? 2 - ^ JV log JJ + (a - 1) JF, (65) 

whose expression matches the one given in [5, 10] when a = \. Since F is a function 
of 9, the above equation is then a non-linear PDE whenever a / |. This is to be 
contrasted with the seemingly much nicer case a = i, whose r.h.s. is independent 

of e. 

A right-hand side independent of 9 can, however, be obtained by taking diver- 
gence of both sides of (65) since F = V log cr and V ■ (JV) = 0. The equation we 
obtain is 

A0 = V-V£. (66) 

This elliptic PDE requires knowledge of 9 at the domain's boundary. Assume that 
we know u,|ox = 9ii Ji = pd^Ui for i = 1,2, and a at the boundary. In this setting, 
we find that 

9\dx = arg^iiVtii + t 12 Vu 2 \ dx ) 

= arg((in<9 t 5i + t 12 d t g2) t + v~ x {twJ\ + t 12 J 2 ) v), 

with u and t = Ju the unit outgoing normal vector and its direct orthogonal vector, 
respectively. 

Once 9 is reconstructed, we know the r.h.s. of (60) and solve for logcr, either by 
integrating (60) along a curve, or by taking the divergence of both sides of (60) and 
solving a Poisson equation provided that cr\gx is known. Note that the inversion 
for 9 and log a by means of the elliptic equations (66) and "divergence of (60)" with 
Dirichlct conditions is unique and Lip schitz- stable in H 2 {X) w.r.t. the data ffy. 
The details are left to the reader. 

We now discuss the compatibility conditions for the gradient equations (60) and 
(65), which admit a solution only if their respective r.h.s. are curl-free. Such 
conditions lead to a better understanding of the range of the measurement operator 
and are necessary to ensure that reconstructions based on ODE integrations do not 
depend on the choice of integration path. 

6.2. Compatibility conditions in two dimensions of space. The compatibil- 
ity conditions for (65) and (60) are that V • (JV0) = and V • (JV logcr) = 0, 
respectively. For a ^ \ , these equations not only provide constraints on the redun- 
dant data, but in fact give us direct information about the unknown coefficients. In 
the two-dimensional case, they allow us to solve algebraically for cos(20), sin(20), 
which in turn characterizes F in terms of the data (and therefore does not require 
the prior resolution of 9). 

Let us first simplify the expression of F as follows: 

F = VlogD + 2($ u Fii +®22V 22 + ($12 + $21)^12) 
= -Vu - V 22 + 2($ u y n + $22^22 + (*i2 + SaOVia) 

= ($11 - $22)(^ll - V22) + ($12 + $2l)(Vl2 + V21), 
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where the matrices $n — $22 and $12 + $21 are reflexion matrices that can be 
expressed in the following manner: 

$11 - $22 = C2U + s 2 JU and $i 2 + $21 = -S2U + c 2 JU, where 



(c 2 ,s 2 ) := (cos(20),sin(20)), U 



1 

-1 



(67) 



(68) 



As a result, we are able to express F in a rather compact way 

F(x, 6) = cos{28)F c {x) + sm(28)JF c {x), where 

f c ( x ) := u(v n - y 22 ) + ju(^ 12 + y 21 ). 

Note the property that dgF = 2JF. Now turning to the compatibility conditions 
proper, equations (65) and (68) arc well-defined only if their curls are zero, which 
gives the following two scalar conditions 

V- (JG) - (a- i)v- (F(x,8(x)) =0 and V ■ ( JF(x, 6(x))) = 0, 

where G := F 12 — - JV log L>. (69) 

Now using the chain rule and V6* = G + [a — i) JF, we have 
V • (JF) = c 2 \7 ■ (JF C ) - s 2 V ■ F c + V6 ■ (JOgF) 

= c 2 V • (JF C ) - s 2 V -F C -2G- (c 2 F c + s 2 JF c ) = fc 2 - gs 2 , 
where we have defined 

f(x) := V ■ [JF C ) — 2F C ■ G and g(x) := V • F c + 2JF C ■ G. (70) 
Similarly, the second compatibility equation can be recast as 

(a - I) (gc 2 + fs 2 ) = V ■ (JG) - 2(a - ^) Vd'- 
If a 7^ h, we thus see that the two compatibility equations imply the system 



/ -9 1 [ c 2 1 = r 

9 f \ [ s 2 J [ft 

which may be inverted as 

cos(26>) = c 2 = 



ft := [a 

2 



V-(JG)-2(a--)|F c | 2 , (71) 



P+9 2 



and sin(2^) = s 2 



/ft 



/ 2 +.9 



2 ' 



Note that this solution makes sense only if the functions /, g, h are such that c 2 +s 2 = 
1, that is, if they satisfy the relation f 2 + g 2 = h 2 . In this case, F may be expressed 
as 



F = V log (J 



( ff F c + / JF C ) . 



/ 2 +.9 2 

The right-hand-side is guaranteed to be curl- free by construction. Inserting (70) 
into the last equation and using (u ■ Jv)v — (u ■ v)Jv — \v\ 2 u with u = G and v = F c , 
we obtain the following explicit reconstruction formula: 

Vloga = (2 k ((V • F C )F C + V ■ {JF C )JF C + 2\F C \ 2 G). (72) 
/ + 9 

When a = |, the equation V ■ (JV0) = depends solely on the data and reads 

V-(J^ Q 2 )-iAlog£ = 0. 
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The other compatibility equation V • ( JF) = is still of the form fc2 + gs 2 = 0, 
with /, g defined above, which by itself only gives us (c2, s 2 ) up to a sign, i.e. 

{c 2l s 2 ) = ±(f+g 2 )- 1 H-gJ). 

The above constraint provides partial answer about 9 that may be used in practical 
reconstructions to mitigate the influence of noise in the data. Reconstructions 
based solely on these algebraic relations, however, seem to be less stable than the 
two approaches based on integration of gradient or Laplace equations. 

6.3. Compatibility conditions in higher dimensions. The two-dimensional 
case is special in that 5*02 (R) is both one dimensional and Abelian. This is not 
the case in higher dimensions, where parameterizations are much more compli- 
cated, even in three dimensions. Similar compatibility conditions arise in the n- 
dimensional case for gradient equations as an application of the Poincare lemma: 
the "curl" (or exterior derivative) of the equation vanishes on both sides. For equa- 
tions (9) or (42), we have that 

dF b = d 2 log a = 0, (73) 

which implies a system of |ri(n — 2) scalar equations (i.e. the cancellation of a 
2-form in dimension n). 

Regarding the systems (35) and (45), their complete integrability is equivalent 
to ensuring that the curvature tensor of the Euclidean connection is identically zero 
when expressed in either frame S or R. Indeed, according to [9, Theorem 1 p30], a 
system of the form 

d k Ej(x)=F ijk (E(x),x), l<i,j,k<n, I:' : { /.';' | ;* ; ; . 
is integrable if and only if the following conditions hold 

diFijk + Fp q id E 9Fijk = dkF lJ i + F pqk d E ^Fiji, 
which is equivalent, after using the chain rule, to 

OidkEj - d k diE{ = di(F m {E(x),x)) - d k {F tjl {E(x),x)) = 0. 

The last equation is nothing but the fact that the curvature (lZ(ei,e k )Ei) ■ ej of the 
Euclidean metric is zero for every quadruple (i,j,k,l), where the curvature tensor 
1Z is defined, for three vector fields X, Y, Z, by 

K(X, Y)Z := VxVyZ - VyV.xZ - V [x ,y]Z. (74) 

1Z is known to be a tensor, i.e. it is linear over smooth functions in all its arguments, 
and thus the above zero curvature conditions are strictly equivalent to the equations 

(1Z(E p ,Eq)Ek) ■ E r = (K{E p e u Ei q e 3 )E k ) ■ E l r ei = E l p E ] q E l J (n(e l ,e 1 )E k ) ■ e, = 0, 

for 1 < p, q, k, r < n, where E is any frame, including S and R. As one can see from 
[21, Prop. 10 and 12 pp 196-197], this highly redundant set of n 4 scalar equations 
is equivalent to ^n(n — 1) non-redundant equations which express the cancellation 
of the sectional curvatures 

{K{Ei, E,j)Ei) ■E j =0, l<i<j<n. (75) 

For each frame, we thus have the two systems of ^n(n — 1) equations (73) and 
(75). We now work with the R frame because its Christoffel symbols have nicer 
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symmetry properties and show that both systems (73) and (75) may be recast as 

n 

M?? : Rp <g> R q = 0, 1 < i < j < n, (76) 

p,q=l 

where the matrices M^ 9 depend only on the data (we could obtain similar equations 
for the S frame). We first rewrite the Christoffel symbols (44) under the form 



^k — u/kl 
ij 

W$ := 2ac F {S ik V l s j - <J y V&) + S u V; k + * y (Vf k + ac F N5 lk ) - S lk (V$ + acpNS^) . 

(77) 

The vector fields Wfj depend only on the data and have the antisymmetry properties 
Wij = —W( k and W{j = 0. We now derive systems of the form (76) for both systems 
(73) and (75) in the case of the R frame. 



rf , = W% ■ R h where 



6.3.1. The condition dF = 0. Starting from equation (42), the equation dF^ = 
reads 

= ^dF" = \d 2 log£> + dm x ■ R k )R\) = d{(V k \ ■ R k )R\). 
Now using identity (26) with vector fields Ri,Rj (i ^ j), we have 
d({V k \ ■ Rk)R\){Ri,Rj) = V B , (V k s 3 ■ R k ) - V Rj {V k \ ■ R k ) - (V& ■ R k )R t ■ [R t ,Rj] 

+ (V k s r R k )((W%-W$)-R p ) = 0. 
Decomposing the first term in the last r.h.s. as follows 

V B , (V£i ■ Rk) = V* \C, ■ Rk + ■ V Ri R k = R k ■ V Ri V^ + (V& ■ R P )(W[ k q ■ R q ), 

and doing simlarly for the second term, we obtain the set of scalar equations 

Rk ■ (V fli V£ - V R] V k s i) + (V£ • R p )(W[ k q ■ R q ) - (V k \ ■ R P ){W™ ■ R q ) 

+ (V q ] ■ R q )((W% - W%) ■R p ) = 0, l<i<j<n. 

This system can be written in the form (76), where the matrices depend only 
on the data. 

6.3.2. The zero curvature conditions. Given the symmetries of the Christoffel sym- 
bols, one can show for an orthonormal frame that the zero sectional curvature 
equations (75) can be recast as 

+ v* r;, = -r'.rl. + - (r^ - r^)4, 1 < i < j < n . (79) 

Using expression (77) of the Christoffel symbols, the first term in the left-hand side 
of (79) may be rewritten as 

v„ r;, = Vr 4 wfj ■ Ri + Wj l j ■ W R /,',. = Ri • V R w% + rjwg • R k 

= R l -V Ri Wf j + {Wf-R p ){Wf r R q ). 
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Proceeding similarly for the second term of the l.h.s. of (79) and plugging expression 
(77) into the r.h.s., we obtain the set of equations for 1 < i < j < n 

Ri ■ (V^j +V Rj W!t) + {Wf ■ R P )(W^ • R q ) + (W™ ■ RJiwg ■ R q ) 

= ■ Rv)(K ■ R «) + ( W u ■ R v)( W n ■ R J ( 8 °) 

+ ((Wtf-W 1 *)-R p )(Wj?-R q ) 

are not being summed over but l,p,q are). This is also a quadratic system of 
the form (76) with different matrices Mt-L 

6.3.3. Discussion. Based on the result of the two-dimensional case, we make the 
following heuristic statements: depending on the value of a, these compatibility 
equations 

(i) either give us compatibility conditions on the data (that do not depend on the 
unknown frame), thus characterizing the range of the measurement operator, 

(ii) or they may allow us to invert algebraically for the cosines and sines of the 
n(n — l)/2 spherical angles that parameterize the SO n (R)-valued R frame. 
This in turn may come at the price of other compatibility conditions that 
only depend on the data. 

It remains an interesting, so far unresolved, question to find an algorithm that 
enforces the compatibility conditions as the system of ODEs is used to ensure that 
the reconstruction does not depend on the choice of integration paths. 

6.4. Remark on the elliptic method. The system of elliptic equations for the 
scalar solutions ui is well-posed for almost all values of a and n. Once the solutions 
Ui are obtained, it remains to solve the equation for F = V log a. The only remain- 
ing compatibility condition is therefore that the latter term indeed be a gradient. 
In a similar manner to what we just saw for the ODE-based method, writing the 
condition dF b = yields \n[n — 1) equations of the type 

M pq : Vu p ® Vu « = °! 1 < i < 3 < 

p,q 

where the matrices depend only on the data Hij . 
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